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ABSTRACT 

We investigate the resonant rotation of circumbinary bodies in planar quasi-circular orbits. Denoting n b and n the orbital mean motion 
of the inner binary and of the circumbinary body, respectively, we show that spin-orbit resonances exist at the frequencies n ± kv/2, 
where v = n b - n, and k is an integer. Moreover, when the libration at natural frequency has the same magnitude as v, the resonances 
overlap and the rotation becomes chaotic. We apply these results to the small satellites in the Pluto-Charon system, and conclude that 
their rotations are likely chaotic. However, the rotation can also be stable and not synchronous for small axial asymmetries. 
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1. Introduction 


2. Spin-orbit coupling 


Circumbinary bodies are objects that orbit around a more mas¬ 
sive binary system. In the solar system, the small satellites of 
the Pluto-Charon system are the best example (e.g., Brozovic 
et al. 2015). Planets orbiting two stars, often called circumbi¬ 
nary planets, have also been reported (e.g., Correia et al. 2005; 
Doyle et al. 2011; Welsh et al. 2012). Most of these bodies are 
close enough to the central binary to undergo tidal dissipation, 
which slowly modifies the rotation rate until it becomes close to 
the mean motion (e.g., MacDonald 1964; Correia et al. 2014). 

For an object that moves in an eccentric orbit around a sin¬ 
gle massive body, the rotation rate can be captured in a half¬ 
integer commensurability with the mean motion, usually called a 
spin-orbit resonance (Colombo 1965; Goldreich & Peale 1966). 
Moreover, for very eccentric orbits or large axial asymmetries, 
the rotational libration width of the individual resonances may 
overlap, and the rotation becomes chaotic (Wisdom et al. 1984; 
Wisdom 1987). However, for nearly circular orbits, all these 
equilibria disappear and the only possibility for the spin is syn¬ 
chronous rotation (e.g., Correia & Laskar 2009). 

When a third body is added to the classic two-body problem, 
the mutual gravitational perturbations introduce additional spin- 
orbit resonances at the perturbing frequency (Goldreich & Peale 
1967; Correia & Robutel 2013). In the case of circumbinary bod¬ 
ies there is, in addition, a permanent misalignment of the long 
inertia axis of the rotating body from the radius vector pointing 
to each inner body (Fig. 1). The resulting torque on the rotating 
body’s figure induces some rotational libration that may give rise 
to some unexpected behaviors for the rotation rate (Showalter & 
Hamilton 2015). In this Letter we investigate all possibilities for 
the rotation of circumbinary bodies, and apply them to the spin 
evolution of the small satellites of the Pluto-Charon system. 


We consider a planar three-body hierarchical system composed 
of a central binary with masses m o and mi, together with an 
external circumbinary companion with mass m (Fig. 1), where 
m <sc mi < mo- For the orbits we use Jacobi canonical coor¬ 
dinates, with rb being the position of mi relative to mo (inner 
orbit), and r the position of m relative to the center of mass of 
mo and mi (outer orbit). The body with mass m has principal 
moments of inertia A < B < C and angular velocity given by a>. 
In Appendix A we provide the full equations of motion for the 
spin. Since tidal dissipation usually damps the obliquity to zero 
(e.g.. Hut 1980; Correia 2009), for simplicity, here we describe 
the motion for a> normal to the orbit. In absence of tides, the 
equation for the rotation angle, 6 , is then given by (Eq. (A. 5)) 
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where a and n are the semi-major axis and the mean motion of 
the outer orbit, respectively, (r,-,/i) are the radial and the angular 
coordinates of r ( . 
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The parameter 6 is the mass ratio of the inner binary, while cr 
is approximately the frequency of small-amplitude rotational li- 
brations in an unperturbed synchronous resonance. 

We can express r 0 3 , sin/o, and cos/o as 


* Fig. 4 and Appendices A and B are available in electronic form at 
http://www.aanda.org 
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Fig. 1. Reference angles and coordinates. The system is coplanar 
and is composed of an inner binary with masses w 0 and nt\, and 
an outer body m that rotates with angular velocity to. 


Each term individually behaves like a pendulum where the ro¬ 
tation can be trapped. Therefore, together with the classic syn¬ 
chronous equilibrium at 6 = n, there are additional possibili¬ 
ties for the spin at the super- and sub-synchronous resonances 
6 = n + kv/2, with k e Z. 

Terms with k + 0 depend on cr 2 , 6, and a 2 , that is, on the 
axial asymmetry of the body, on the mass ratio of the central bi¬ 
nary, and on the semi-major axis ratio. Thus, for circumbinary 
bodies far from the inner binary (a <s 1), the amplitude of these 
higher order terms decreases very quickly to zero. For <5 « 0 (i.e., 
m i <s mo) we have P 2 - Pi = 0, thus only the synchronous reso¬ 
nance subsists, as in the classic circular two-body problem (e.g., 
Goldreich & Peale 1966). The effect of the binary mass ratio is 
maximized for 6 - 1/2 (i.e., mo = m{), but only for the second 
order resonances (terms in p 2 ), since third order resonances also 
vanish for equal masses (p 3 = 0). Interestingly, the size of the 
sub-synchronous resonance with k < 0 is always larger than the 
symmetrical super-synchronous resonance with k > 0. 


and 
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Similar expressions can be found for r/ 3 , sin/), and cos/), if 
we replace 6 by (6 - 1). Therefore, since ||rj,|| < ||r||, we can 
develop equation (1) in power series of r^/r. In addition, for the 
hierarchical problem the orbits are approximately ellipses. Thus, 
r and r* can also be expanded in terms of Hansen coefficients 
(e.g., Laskar & Boue 2010). We then rewrite equation (1) as 

cr 2 v— 1 

0 = —y Zj eb ’ ^ Sln ~ P M - c i M ^ > ( 6 ) 

lp,q 


where e is the eccentricity; M is the mean anomaly; Pi, IKq °c a 1 , 
with a = cib/cr, and p, q are half-integers. Spin-orbit resonances 
occur whenever 6 — pn + cpih. In Appendix B we provide the 
complete expression for equation (6) to the first order in the ec¬ 
centricities and to the second order in a. 


3. Near circular orbits 

For near circular orbits (e * 0 and e/, ~ 0), we can neglect both 
eccentricities, thus r « a and /-/, ~ a*. Retaining terms in a 3 in 
equation (6) gives 
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where y — 8 - M, tf> = Mf, - M, 

p 2 =5(l-6)a 2 , and p 3 = 6(1 - 6)(1 - 26) a 3 . (8) 

We see that we have several islands of rotational libration, 
y — Q - n — 0, +v, ±v/2, ±3v/2 , with v = 0 = rib - n . (9) 


3. 1. Chaotic motion 

When the libration amplitudes of some individual resonant is¬ 
lands overlap, the rotational motion can be chaotic (Chirikov 
1979; Morbidelli 2002). This means that the rotation exhibits 
random variations in short periods of time. This phenomenon 
was first described for the rotation of Hyperion, which is chaotic 
because its orbit is eccentric (Wisdom et al. 1984). 

In the circular approximation, each individual resonance is 
placed at y — kv/2, with libration width cr\[Pk, where //j = 
YjiPuo-k)l 2 ,ki 2 - Overlap between two resonances with k\ < k 2 
then occurs whenever 

k\v/2 + cr > k 2 v/2 - cr y[Pk 2 . (10) 

The synchronous resonance (k = 0) has the largest width, so 
overlap is more likely for this resonance. The third order reso¬ 
nance with k - -1 is the nearest resonance with larger ampli¬ 
tude, thus chaos sets in for 


cr 



1 + 7 P 2 



2 («6 - n) . 
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Rotation in the chaotic zone is also attitude unstable 
(Wisdom et al. 1984; Wisdom 1987). A non-zero obliquity intro¬ 
duces additional resonant terms to the rotation (Eq. (7)), increas¬ 
ing the chances of overlapping (e.g., Correia & Faskar 2010). 
Therefore, in order to correctly account for the chaotic behavior, 
we need to integrate the full equations of motion (Appendix A). 


3.2. Global dynamics 

The general problem of the spin-orbit coupling for near circular 
circumbinary bodies can be reduced to the analysis of y/n, which 
depends on the libration frequency cr/n, the mass ratio 6, and the 
semi-major axis ratio a (Eq. (7)). However, since <5 e [0,1/2] 
and it is a well-known parameter for a specific binary system, 
the global dynamics is approximately controlled by the only two 
free parameters: a (or nb/n ) and cr/n (or (B - A)/C). 

Rescaling the time with n = 1, we can perform a stability 
analysis of y in the plane (a, cr) to quickly identify the rotational 
regime for any system of quasi-circular circumbinay bodies that 
is near the synchronous equilibrium. If isolated, the half-width of 
the synchronous resonant island in the direction of y is equal to 
cr. When perturbed, chaos is expected for cr ~ v (Eq. (11)). Thus, 


2 

















Alexandre C. M. Correia et al.: Spin-orbit coupling for circumbinary bodies 



o 0.2 0.4 0.6 0.8 1 

Fig. 2. Stability analysis of the rotation rate close to the synchro¬ 
nization for 5 = 0.1, a e [0.25 : 0.5], and (B - A)/C e [0 : 0.9]. 
We adopt three different initial obliquity values i// = 1°, 30°, and 
60°. The color index indicates the proportion of chaotic orbits 
inside the studied domain: from dark blue for fully regular to red 
for entirely chaotic. The dashed line (Eq. (11)) roughly delimits 
the zone where chaotic motion is possible. The dots identify two 
representative cases that we detail in Figure 3. 


for a given (a, u ) pair, we fix the initial value of y at the libration 
center of the synchronous resonance, and select 400 equispaced 
initial values of y e [-v : v]. The corresponding solutions are in¬ 
tegrated using the equations for the three-body problem together 
with equations (A. 1) and (A.2) for the rotational motion (without 
tidal dissipation). The dynamical nature (stable/unstable) is de¬ 
duced from frequency analysis (Laskar 1990, 1993), which gives 
the fraction of chaotic trajectories. 

In Figure 2 we show the general results for 6 - 0.1 using 
three different initial obliquity values i// = 1°, 30°, and 60°. All 
the remaining variables are initially set to zero. The color index 


indicates the proportion of chaotic orbits inside the studied do¬ 
main: from dark blue for fully regular to red for entirely chaotic. 
Depending on the value of a and cr, the rotation can present 
different behaviors, ranging from non-synchronous equilibria to 
chaotic motion. Figure 2 remains almost unchanged for <5 > 0.1 
since the global dynamics is not very sensitive to 6. The only 
exception would be for 6 <s 0.1 since in that case the system 
would behave like the classic circular two-body problem where 
the rotation can only be synchronous. 

A common way of visualizing and understanding the differ¬ 
ent regimes is to use frequency map analysis (Faskar 1993). We 
selected two representative pairs (a, cr) and plotted the corre¬ 
sponding cross-section maps for the three initial obliquities in 
Figure 3. We are able to identify the stable and the chaotic re¬ 
gions very easily. We observe that stable islands can exist in the 
middle of chaos, implying a large diversity of possibilities for the 
final rotation. In the chaotic regions for the rotation, the obliquity 
is also chaotic and can vary between 0° and 180°. For high obliq¬ 
uities, the size of the stable synchronous island shrinks and the 
chaotic zone is extended. Therefore, depending on the obliquity, 
the rotation alternates between more or less chaotic. In times of 
small obliquity, the rotation can be trapped in a stable resonant 
island, preventing the obliquity from increasing again, and the 
spin can become stable (see Fig. 4). 

For small values of a and/or cr, the individual resonances 
(Eq. (7)), associated with the plateaus in Figure 3, are well sepa¬ 
rated. Thus, for dissipative systems, the rotation can be captured 
in individual spin-orbit resonances and stay there. For rotation 
rates increasing from lower values, the sub-synchronous reso¬ 
nance with 9 = n- v = 2n — rib has the largest amplitude, so 
this is the most likely spin-orbit resonance to occur. For rotation 
rates decreasing from higher values, the super-synchronous res¬ 
onance with 6 = n + v = iib becomes the most likely possibility. 
This case is very interesting as it corresponds to a synchronous 
rotation with the inner orbit period. However, since the ampli¬ 
tude of this resonance is small, the rotation can easily escape 
it and subsequently evolve into the classic synchronous rotation 
with the outer orbit (9 = n). 

The amplitude of the resonant terms increases with a and cr. 
At some point, the individual libration islands merge and chaotic 
motion can be expected (see section 3.1). The transition between 
the two regimes is roughly given by the dashed curved, obtained 
with expression (11). Near the transition regime, chaotic rotation 
can be observed around the separatrix of the synchronous reso¬ 
nance, but the motion is still regular nearer the center (Fig. 3, 
left). Therefore, when the rotation of a body is evolving by tidal 
effect, the rotation becomes chaotic when approaching the syn¬ 
chronous resonance. However, depending on the strength and on 
the geometry of the tidal torque, it is still possible to find a path¬ 
way to the synchronous rotation. The wandering in the chaotic 
region may also provide a path into another stable super- or sub- 
synchronous resonance (see Fig. 4). 

For a > 0.4 and cr ~ 1, there is a large overlap between sev¬ 
eral resonant terms, in particular for the negative ones (k < 0). 
As a consequence, there is a large chaotic zone for the spin 
(Fig. 3, bottom). Some of the individual resonances may still 
subsist, including the synchronous one, but they have small sta¬ 
ble widths and they can only be reached at periods of small 
obliquity. The chaotic motion is maximized for cr « 1, which 
corresponds to (B - A)/C ~ 0.35, because secondary reso¬ 
nances generate strong instabilities inside the synchronous is¬ 
land (Robutel et al. 2012). Finally, for a ~ 0.49 the system 
reaches the 3/1 mean motion resonance, which introduces addi¬ 
tional forcing to the rotation, and all trajectories become chaotic. 
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Fig. 3. Cross sections representing the spin dynamics of Hydra (top) and Nix (bottom) at different initial obliquities ip = 1°, 30°, 
and 60°. The dots show the main frequency of y(t), denoted by for different initial values of y(0). It is obtained by numerical 
integration of the three-body problem together with equations (A.l) and (A.2) for the rotational motion (without dissipation). 6 is 
initially set equal to / so the sections go through the middle of the synchronous resonance (j = 0). All the remaining variables 
are initially set to zero. The dots are blue for stable rotation (libration or circulation), and red for chaotic motion. The distinction 
between stable and unstable trajectories is based on the estimate of the second derivative of with respect to the initial condition 
y(0) (see Laskar 1993). The plateaus correspond to the resonance crossings (£/v = k/2). The red dots randomly distributed between 
the plateaus indicate the overlapping of the associated spin-orbit resonances, which generates chaotic motion. 


4. Application to the Pluto-Charon system 

In 1978, regular series of observations of Pluto showed that 
the images were consistently elongated revealing the presence 
of Charon (Christy & Harrington 1978). The orbital parameters 
hence determined have shown that the two bodies evolved in an 
almost circular orbit with a 6.4-day period (e.g., Tholen & Buie 
1997). Charon has an important fraction of the total mass (about 
11%), and therefore the system is considered a binary planet 
rather than a planet and a moon. Later, it was found that four 
tiny satellites move around the barycenter of the Pluto-Charon 
system (Weaver et al. 2006; Brozovic et al. 2015), also in nearly 
circular and coplanar orbits (see Table 1). 

Pluto’s brightness also varies with a period of 6.4 days (e.g., 
Tholen & Tedesco 1994). Since Charon is too dim to account for 
the amplitude of the variation, this period has been identified as 
the rotation of Pluto. Therefore, at present, the spin of Pluto is 
synchronous with the orbit of Charon, a configuration acquired 
from the action of tidal torques raised on Pluto by Charon (e.g., 
Farinella et al. 1979; Cheng et al. 2014). Tidal torques raised 
on Charon and on the remaining satellites by Pluto are even 
stronger, so all these bodies are presumably also tidally evolved. 


Table 1. Mean orbital and physical parameters for Pluto’s moons 
(Brozovic et al. 2015). The parameter mo is the mass of Pluto. 


Par am. 

Charon 

Styx 

Nix 

Kerberos 

Hydra 

P [day] 

6.3872 

20.1617 

24.8548 

32.1679 

38.2021 

n b /n .. 

1 

3.1565 

3.8913 

5.0363 

5.9810 

a .... 

1 

0.46203 

0.40246 

0.33932 

0.30278 

e . 

0.00005 

0.00001 

0.00000 

0.00000 

0.00554 

m/mg . 

0.12176 

oio - 7 

3.4 10- 6 

1.3 10~ 6 

3.7 10~ 6 

R [km] 

604 

4-14 

23-70 

7-22 

29-86 


Showalter & Hamilton (2015) have measured the brightness 
variations of the small satellites. They concluded that Hydra and 
Nix show no obvious pattern, suggesting that their rotation is 
chaotic. Assuming that these two satellites are uniform triaxial 
ellipsoids, they additionally estimate their figure semi-axis ra¬ 
tios, which allows us to compute B/C ~ 0.92 and A/C w 0.44 
for Hydra, and B/C ~ 0.94 and A/C ~ 0.31 for Nix. 

In Figure 2 we already show the global dynamics for the 
Pluto-Charon system, for which 6 = 0.1085 « 0.1 (Brozovic 
et al. 2015). Therefore, we plot vertical lines at the a values 
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corresponding to the small satellites in the system (Table 1). 
Moreover, since we have an estimation for the axial asymme¬ 
tries of Hydra and Nix, we also show a dot for their ( B - A)/C 
values in Figure 2. These two satellites are in a different dynam¬ 
ical regime: the rotation of Hydra can present a stable spin-orbit 
coupling (Fig. 3, top), while the rotation of Nix is most likely 
chaotic (Fig. 3, bottom), as discussed in section 3.2. The spin dy¬ 
namics of Kerberos is probably similar to that of Hydra, while 
the rotation of Styx can be even more chaotic than that of Nix, 
depending on their (B - A)/C values. 

To test the reliability of the dynamical picture described 
in the previous section, we now add the effect from tides to 
our model, so that we can follow the long-term evolution of 
the spin of these satellites. For simplicity, we adopt a constant 
time-lag linear model, whose contribution to the spin is given 
by expression (A. 13). In Figure 4 we show some examples for 
the final evolution of Hydra’s spin with tidal dissipation, start¬ 
ing with slightly different initial values of 0. We adopt an ini¬ 
tial retrograde rotation of 4.4 days (y/v = -2.3), and 30° for 
the initial obliquity, to force the rotation to cross the large am¬ 
plitude sub-synchronous resonances. An initial retrograde rota¬ 
tion for such small bodies is as likely as a prograde one (e.g., 
Dones & Tremaine 1993). For the uncertain parameters we used 
R = 45 km, k 2 /Q = 1(T 4 , and C/(mR 2 ) = 0.4. 

In one example (Fig. 4a) the rotation of Hydra is trapped in 
the spin-orbit resonance with y = -3v/2, while in the others the 
rotation reaches the chaotic zone. However, after some wander¬ 
ing in this zone some simulations are able to find a path into a 
stable spin-orbit resonance with k > 0 (Fig. 4b) or with k < 0 
(Fig. 4c). Interestingly, while for resonances with k > 0 the spin 
axis stabilizes near zero degrees, for resonances with k < 0, the 
spin axis stabilizes with a high obliquity value. The last example 
(Fig. 4d) remained chaotic for the length of the integration and 
therefore can represent the observed present state. 

Over 10 Myr, most of our simulations remained chaotic, 
but several captures in stable non-synchronous resonances also 
occurred. The final scenario depends on the initial conditions, 
and also on the tidal model. A constant-Q model would pre¬ 
vent any capture in resonance (Goldreich & Peale 1966), while a 
visco-elastic model would increase the chances of capture (e.g., 
Makarov 2012; Correia et al. 2014). Strange attractors can also 
exist in the chaotic zone, which may prevent the spin from sta¬ 
bilizing (e.g., Batygin & Morbidelli 2011). However, the global 
dynamics described in section 3.2 is very robust and reliable, 
since it does not depend on dissipative forces (Fig. 2). 

5. Discussion 

This work was motivated by the recent observations on the ro¬ 
tation of Hydra and Nix (Showalter & Hamilton 2015), which 
appear to be chaotic. Our model confirms that chaotic rotation 
is a likely scenario for both satellites, but stable spin-orbit cou¬ 
pling could also be possible, in particular for Hydra. This model 
assumes a hierarchical three-body system with coplanar orbits. 
However, we also integrated the spins of Hydra and Nix using 
the ephemerides for the full system provided by Brozovic et al. 
(2015) and SPICE routines (Acton 1996), and we have found no 
evidence of any substantial differences. 

In our study, we conclude that stable spin-orbit coupling is 
also a plausible scenario for near circular circumbinary bodies 
with small a and/or small <x (Fig. 2). Equilibrium rotation occurs 
for 6 = n+kv/2, with v = H),-n and k e 2. The largest amplitude 
non-synchronous resonance corresponds to the sub-synchronous 
resonance at n - v = 2 n - . Bodies captured in this resonance 


present retrograde rotation if rib/n > 2. This condition is verified 
for the small satellites of the Pluto-Charon system (Table 1), and 
likely for any circumbinary system, since large orbital instabil¬ 
ities are expected for period ratios below the 2/1 mean motion 
resonance. Therefore, bodies trapped in a non-synchronous res¬ 
onance are also likely to present retrograde rotation. 

Lately, many planets have been detected around binary stars 
(e.g., Welsh et al. 2012). So far, all these planets are gaseous 
giants, for which the axial asymmetry is very low, for instance 
( B - A)/C ~ 10 7 for Jupiter (Jacobson 2001). Therefore, al¬ 
though many of these planets are close enough to their stars to 
undergo tidal dissipation, spin-orbit coupling is very unlikely. 
However, for smaller mass Earth-like circumbinary planets these 
states are possible since (B - A)/C ~ 10 -5 (e.g., Yoder 1995). 
In particular, non-synchronous rotation is possible, which is an 
important point to take into account in future habitability stud¬ 
ies (e.g., Selsis et al. 2007). Unlike the Pluto-Charon system, 
circumbinary exoplanets usually present eccentric orbits. As a 
consequence, the number of spin-orbit resonances drastically in¬ 
creases (Eq. (B.4)). Overlap of the different contributions is then 
easier, so chaotic rotation can also be more likely in this case. 
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Fig. 4. Possible evolutions for the spin of Hydra. We show the rotation y/v (in red) and the cosine of the obliquity (in blue). We 
numerically integrated the equations of the three-body problem together with equations (A.l), (A.2), and (A. 13) for the rotational 
motion. We adopt an initial retrograde rotation of 4.4 days with different initial values of 6 , and an initial obliquity of 30°. 
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Appendix A: Equations of motion 

We let ( a,b,c ) be a non-inertial frame attached to the body’s 
principal inertial axes, with inertia tensor I = diag(A, B, C). In 
this frame, 10 - {a> a ,u)b,<jj c ) and the angular momentum L = 
(Ato a , Bub, Ca> c ). The equations of motion are given by 

— Tq-\-T\ , (A. 1) 


Appendix B: Series with eccentricity 

Equation (1) for the rotational motion can be expanded in power 
series of the eccentricities and semi-major axis ratio a as given 
by the general expression (6) 



where T) is the gravitational torque on the body’s figure. For 
second order interactions, we have (e.g., Goldstein 1950) 

T, = [(^ -A)yirjXb + (C-A) n r, x c] , (A.2) 


where r, has coordinates (v,,y,, z,) in the body’s frame. Thus, 
projecting equation (A.l) over each axis (a, b , c) gives 
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(A.3) 

(A.4) 

(A.5) 


To solve these equations, a set of generalized coordinates to 
specify the orientation of the satellite must be chosen. We adopt 
the modified Euler angles ( 6 , i p, if/) as defined in Wisdom et al. 
(1984). Starting the c axis coincident with the normal to the or¬ 
bital plane, and a along an inertial direction, we first rotate the 
body about the c axis by an angle 0 , then we rotate about a by an 
angle p>, and finally we rotate about the b by an angle if/. Then, 


xjrj = cos (6 - fi) cos if/ - sinffl - fi) sin if/ sin tp , (A.6) 

yjn = - sin(0 - fi) cos tp , (A.7) 

Zi/ri = cos(0 - f) sin i f + sin(0 - f) cost//simp , (A.8) 

and 

8 - (to c cos if/ - cj a sin if/)/ cos tp , (A.9) 

ip — a> c sin f + <o a cos f , (A. 10) 

f — a>b - (oj c cos f - oj a sin f) tan tp . (A.ll) 


The full equations of motion for the spin are thus described 
by the set of variables (u> a , a>b, co c , 6, tp, if/), whose derivatives are 
given by equations (A.3)-(A.5) together with (A.9)-(A.l 1). The 
angle e between the c axis and the normal to the orbit, usually 
called obliquity, can be obtained from 


with 

T = Pi,p,q(e, e b , 6, a) sin 2(0 - pM - qM b ) . (B.2) 

When we truncate the series to the first order in the eccentricities 
and to the second order in a, we get 


T = (l + -p 2 ) sin(20 - 2M) 
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(B.3) 


cos s = cos tp cos if/ . 


(A. 12) 


For the dissipative tidal torque, we adopt a constant time-lag 
linear model, whose contribution to the spin (Eq. (A.l)) is given 
by (e.g., Mignard 1979) 

2 

Td = 3k 2 GR 5 At ^ ~y [(/•,■ ■ cj) r, - r 2 o/ + r, x r,] , (A. 13) 

1=0,1 r i 


where is the Love number. At is the time lag, Q 1 = nAt is the 
dissipation factor, and R is the radius of the rotating body. Also 
note that r, is the derivative of r, in an inertial reference frame. 
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